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Abstract. We present a Monte-Carlo calculation of the microcanonical ensemble of the of the ideal hadron- 
resonance gas including all known states up to a mass of about 1.8 GeV and full quantum statistics. 
The microcanonical average multiplicities of the various hadron species are found to converge to the 
canonical ones for moderately low values of the total energy, around 8 GeV, thus bearing out previous 
analyses of hadronic multiplicities in the canonical ensemble. The main numerical computing method is 
an importance sampling Monte-Carlo algorithm using the product of Poisson distributions to generate 
multi-hadronic channels. It is shown that the use of this multi-Poisson distribution allows an efficient and 
fast computation of averages, which can be further improved in the limit of very large clusters. We have 
also studied the fitness of a previously proposed computing method, based on the Metropolis Monte-Carlo 
algorithm, for event generation in the statistical hadronization model. We find that the use of the multi- 
Poisson distribution as proposal matrix dramatically improves the computation performance. However, 
due to the correlation of subsequent samples, this method proves to be generally less robust and effective 
than the importance sampling method. 
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1 Introduction 

The basic assumption of the statistical hadronization model 
(SHM) is that, as a consequence of a dynamical QCD- 
driven process, the final stage of a high energy collision 
gives rise to the formation a set of colourless massive 
extended objects, defined as clusters (or fireballs). They 
are assumed to produce hadrons in a purely statistical 
manner, namely all multihadronic states within the clus- 
ter volume and compatible with its quantum numbers are 
equally likely. The set of equiprobable states with fixed 
four-momentum and internal charges (abelian or not) is 
what is usually called the microcanonical ensemble . In- 
teractions between stable hadrons are mostly taken into 
account by the inclusion of all resonances as independent 
states 1,, which is the reason of the usual expression, in 
the SHM, of ideal hadron-resonance gas 0. 

The special motivations for a detailed study of the mi- 
crocanonical ensemble of the ideal hadron-resonance gas 
have been extensively discussed in the first paper {3] . They 
can be shortly summarized as: 

1. The need of a tool to hadronize final-state clusters 
in particle and heavy ion collisions in the statistical 



1 The full microcanonical ensemble should in principle in- 
clude angular momentum and parity conservation. In this 
work, as well as in previous literature, angular momentum and 
parity are disregarded and the resulting set of states are still 
defined as microcanonical ensemble. 



hadronization model and test observables which can- 
not be calculated analytically. 
2. The assessment of the validity of previous calculations 
in the canonical ensemble, especially in the analysis of 
average multiplicities in high energy elementary colli- 
sions imnirj] . 

The first point is quite clear: if we could compute 
microcanonical averages numerically in a practical way, 
we would be able to make predictions on many observ- 
ables within the basic framework of the statistical model 
without invoking further assumptions or approximations 
which are needed to obtain analytical expressions. Furter- 
more, the availability of a Monte-Carlo integration algo- 
rithm opens the way to event generators with hadroniza- 
tion stage modelled by SHM. The second point is some- 
how related to the first. Indeed, in previous work within 
SHM, specific assumptions were invoked in order to al- 
low the use of the canonical ensemble, which is far easier 
to handle with regard to both analytical and numerical 
calculations. It is then necessary to verify whether it was 
sound to calculate multiplicities in the canonical ensemble 
by comparing them to those in the microcanonical ensem- 
ble with the same values of total mass and volume. 

The microcanonical formalism for the hadron gas and 
the relation between microcanonical and canonical ensem- 
bles have been developed in the first paper .'i . In this 
paper we will focus on the numerical analysis; we will dis- 
cuss the main computational issues and we will describe 
two algorithms to sample the microcanonical hadron gas 
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phase space. As has been mentioned above, we will in- 
spect the differences between microcanonical and canon- 
ical averages. In this first study, we will only deal with 
observables pertaining to particle multiplicities, leaving 
the analysis of momentum spectra to further works. 

The microcanonical ensemble of the hadron gas has 
been investigated numerically by K. Werner and J. Aiche- 
lin [7] by using a Monte-Carlo method based on Metropo- 
lis algorithm. Results on specific observables have recently 
been published [Hj in a hadron gas model including fun- 
damental multiplets (two meson nonets and baryon octet 
plus decaplet). In comparison with these previous works, 
we will show calculations for the full hadron gas includ- 
ing all known species (more than 250) up to a mass of 
about 1.8 GeV and, particularly, we will introduce a new 
updating rule for Metropolis algorithm leading to a dra- 
matic improvement of its performance in terms of comput- 
ing time. Moreover, we propose a different Monte-Carlo 
computing method, based on the importance sampling of 
multi-hadronic channel space, which proves to be more 
effective than Metropolis algorithm for the calculation of 
averages. 

The paper is organized as follows: the basic micro- 
canonical formalism is summarized in Sect. 2; the numeri- 
cal method to compute the phase space volumes for a given 
multi-hadronic channel is discussed in Sect. 3; in Sect. 4 
we describe the importance sampling method which is well 
suited to compute averages in the microcanonical ensem- 
ble of the ideal hadron-resonance gas; in Sect. 5 we show 
the comparison between microcanonical and canonical av- 
erages; in Sect. 6 and 7 the Metropolis algorithm is studied 
in detail; finally conclusions are summarized in Sect. 8. 



2 Microcanonical partition function 

In principle, any average on a given statistical mechan- 
ics ensemble could be calculated from the partition func- 
tion. The microcanonical partition function of the hadron 
gas is best defined 3 as the sum over all multihadronic 
states localized within the cluster \h\r) constrained with 
four-momentum and abelian (i.e. additive) charge conser- 
vation: 



n = J2(hv\S\P - Po P )S QM Jhy) 



(1) 



where Q = (Qi, . . . , Qu) is a vector of M integer abelian 
charges (electric, baryon number, strangeness etc.), P the 
four- momentum of the cluster and P op , Q op the relevant 
operators. Provided that relativistic quantum field effects 
are neglected and the volume of the cluster is large enough 
to allow the approximation of finite-volume Fourier inte- 
grals with Dirac deltas, it can be proved [Hj that Q can be 
written as a multiple integral: 
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where qj is the vector of the abelian charges for the j th 
hadron species, Jj its spin, V the volume of the cluster; the 
upper sign applies to fermions, the lower to bosons. The in- 
tegral J2J) is more easily calculable in the rest frame of the 
cluster where P = (M, 0). Unfortunately, an analytical 
solution with no charge constraint (the so-called grand- 
microcanonical partition function) is known only in two 
limiting cases: non-relativistic and ultra-relativistic (i.e. 
with all particle masses set to zero). The full relativistic 
case has been attacked with several kinds of expansions 
9 but none of them proved to be fully satisfactory as 
the achieved accuracy in the estimation of different kinds 
of averages could vary from some percent to a factor 10. 
Therefore, a numerical integration of Eq. (J2J is needed. 
The most suitable method is to decompose f2 into the 
sum of the phase space volumes with fixed particle multi- 
plicities for each species: 



Q — 2J ^{Nj}' 



(3) 



{Nj} being a vector of K integer numbers {N\, . . . , Nk), 
i.e. the multiplicities of all of the K hadronic species. The 
set of integers {Nj} is also defined as a channel because it 
characterizes a specific decay channel of the cluster. The 
general phase space volume ^{n,} for the channel {Nj} 
obtained from can be written as a cluster decomposition. 
Let j be the integer index running over all hadron species, 
and {h nj } a partition (relevant to the species j) of Nj in 
the multiplicity representation, i.e. a set of integers such 
that Nj = J2n-=i njh nj ; let Hj = E„/=i h n } and let c tj 
be the cyclic permutations of the first nj . integers deter- 
mined by the partition, with ^2i =i n ij ~ Nj. Provided 
that relativistic quantum field effects are neglected (which 
is possible if the cluster linear size is greater than pion 
Compton wavelength, 1.4 fm), the phase space volume for 
fixed multiplicities reads J3J: 

n {Nj} = J d 3 Pl ...d 3 Pw s 4 (p-p f ) 

>nE ( ^r:'" fK «> 



{'•... \ II, .", / '» 



where Pt is the sum of all particle four-momenta. The 
factors F ni in the equation above are Fourier integrals 
over the cluster region with proper volume V: 



f % = n ?2^) 



d x e 1 i i (5) 



and the p's are the particle momenta. For sufficiently large 
V, the integral in Eq. JSJ tends to a product of Dirac delta 
distribution and, if we use this limit in Eq. (T4J , we arrive 
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at this expression of Q^^: 
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where, for a set of partitions {h ni }, . . . , {fon K } for each 
of the hadron species, the four-momenta p[. are those of 
lumps of particles of the same species j (Hj in number) 
with mass rijirij and spin Jj. 

For sufficiently large volumes, the leading term in Eq. 
10, henceforth defined as f2^ N y, is that with the maximal 
power of V, i.e. with Hj — Nj Vj. This term corresponds 
to the partitions {h nj } — (Nj,0,...) Vj, namely to one 
particle per lump, and reads: 



{Nj} 



n 



V N ^(2J 1 + 1) 
(2n) 3N JNj\ 
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where N = 52 • AT,- . This is indeed the phase space vol- 
ume in the classical Boltzmann statistics. Therefore, the 
terms beyond the leading one Q in the expansions J2J 
and @) account for the quantum statistics effects. The 
very same expression Q holds as the leading term of the 
more general cluster decomposition Q. In fact, the cyclic 
permutations cu corresponding to the partition {h nj } = 
(Nj, 0, . . .) are identities, thus implying, according to Eq. 



Ni 



n ^ - 



V 
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In the present work, we have used the approximated ex- 
pression JfjJ of Eq. (TJJ to evaluate the phase space volume 
for fixed multiplicities. For the considered cluster masses 
and volumes (see later on) this approximation is satisfac- 
tory for most observables, taking into account that the 
leading term Q of the cluster decomposition is the same 
in both Eq. and © and that subleading terms give at 
most a 10% correction to the leading term. 

A nice feature of the cluster decomposition ijfjj is that 
every term of the expansion is an integral just like the 
classical phase space volume Q with lumps replacing par- 
ticles. Specifically, the Eq. © can be rewritten as: 



Q 



{Nj} 



{h ni },...,{h„ K } 



n 



(Tl 



iNj+H,- 



TjNj 
lln 3 =l 



{Hj} 



(9) 



Note that the factor fln^i h*h ' nas been absorbed in 
flfff .x as it takes into account the identity of the lumps. 
This form of the cluster decomposition shows that in ac- 
tual numerical calculations all of the terms can be com- 
puted with the same routine. 



As has been mentioned in the introduction, in this pa- 
per we are mainly interested in the calculation of quan- 
tities relevant to particle multiplicities and not to their 
momenta, namely their kinematical state. The average of 
an observable O depending on particle multiplicities in the 
microcanonical ensemble can then be written as: 

(0) = Ew} ^})^-}^^ (10) 

Altogether, what we need to calculate in order to evaluate 
an average (|10|l are integrals like ||BJ and Q . The descrip- 
tion of a suitable numerical technique to do that is the 
subject of next section. 



3 Numerical calculation of the phase space 
volume 

In order to calculate efficiently and quickly the phase space 
volume for fixed multiplicities ^{n 3 -} i n Eq. ©, we have 
adopted a Monte-Carlo integration method proposed by 
Cerulus and Hagedorn in the 60's [9,10 and later em- 
ployed by Werner and Aichelin [7J- The method is de- 
scribed in detail in ref. [7] . As we made only slight modifi- 
cations, we just sketch it here; a more detailed description 
is given in Appendix A. As has already been mentioned, 
every term in the cluster decomposition JfjJ is an integral 
of the kind (JJJ and can be calculated by the same numer- 
ical method. 

The calculation is carried out in the cluster's rest frame, 
where P — (M, 0) and V is the proper volume. The inte- 
gral is first written as the product: 



(8) ^{Nj} 



(2ir) 3N 11 A~! 



$(M,mx,...,m N ) 
(11) 

where T = M — 52i=i TO * ^ s ^ ne available total kinetic 
energy for the N particles (or lumps), and <P is an adi- 
mensional kinematic integral: 



<2> 



l r N 

J d 3 Pl . . . dV <5 4 (P - (12) 



After a sequence of variable changes, the function <P is 
rewritten as an integral of an adimensional function T of 
N — 1 variables r,; 6 [0, 1] (sec Appendix A): 

*= / d ri ... [ dr N -xr(n,...,r N -i) (13) 
Jo Jo 

which can be estimated through Monte-Carlo integration 
as: 

N s 



$ = —Yr(r{ 



(*) Jk) 



(14) 



k=l 



(k) 

where r\ are random numbers uniformcly distributed in 
the interval [0, 1] and Ns is the number of samples. 
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Fig. 1. Distribution of 1?{jv.} values computed with a nu- 
merical Monte-Carlo integration based on 1000 random sam- 
ples, along with gaussian fits, for a cluster of 4 GeV mass 
and 0.4 GeV/fm 3 energy density. Three channels are shown: 
(a) 7r°,7r°,7r + ,7r~; (b) 7r°, 7r°, ty + , ir + , ir~ , tt~; (c) 77,77,77,77. Top 
plots refer to the full quantum statistics calculation while those 
below to the classical statistics approximation. 



Fig. 2. CPU time needed to calculate the phase space vol- 
ume of a channel with N neutral pions in a cluster with mass 
N(m v a +500 MeV) and energy density 0.4 GeV/fm 3 . The full 
dots refer to classical statistics and the hollow dots to full quan- 
tum statistics; lines are drawn to guide the eye. The CPU time 
has been normalized to a Pentium IV processor with 2 GHz 
clock rate. 



In order to calculate the full cluster decomposition of 
the phase space volume, in Eq. ©, the above calculation 
is repeated for all of the terms. To reduce the number of 
calls to the random number generation subroutine, one 
can take advantage of rewriting Eq. © as: 



a 



{h ni },...,{h nK } 

yHrp3H-± 



n 



( T l)JV J +g J ( 2 j j 

11 r 



(27T) 



3H 



#(AT,mi,...,Mh) 



(15) 



where fii are meant to be the masses of the lumps defined 
by the partitions {h nj } and H = VV Hj. Since H < N, 
the first H — 1 random numbers, out of N — 1 extracted, 
can be used to estimate the <P integral, according to Eq. 
<jTf|) . for each term in Eq. lfT5)l. 

For Ns = 1000 random samples, the accuracy of the 
Monte-Carlo integration is satisfactory for all channels 
and turns out to be of the order of some percent. This 
may not be sufficient if one is interested in accurate cal- 
culations of single channel rates themselves, but is very 
good for the estimate of other, more inclusive, averages 
such as mean multiplicities, multiplicity distributions etc., 
in which the errors on different channels add incoherently 
(see Eq. (|10f) ') and the relative statistical error is therefore 
greatly reduced with respect to that on the single term 
/2|jVj}- In fig. we show the gaussian distributions ob- 



tained by repeating 10000 times the calculation of ^{jv^} 
for three different channels in a cluster with 4 GeV mass 
and 0.4 GeV/fm 3 energy density. 2 The plots on the top 
row refer to the f2{N } value calculated in full quantum 
statistics, according to Eq. ©, and those on the bottom 
row in classical statistics, i.e. retaining only the leading 
term {7J . It can be seen that the central values in full quan- 
tum statistics are higher than their classical approxima- 
tions, as expected for channels involving identical bosons. 
The calculation in full quantum statistics with the exact 
expression J2J for finite volume yields results very close to 
those obtained with the approximated one © . The spread 
in value of the statistical error (see resolutions quoted in 
fig. 2]) is owing to the different multiplicities and masses 
of the particles in the channels. However, this spread is 
fairly small and stays well below an order of magnitude 
so that, in fact, a fixed number of Monte-Carlo samples 
N$ is appropriate to calculate most fi^jy^'s with a given 
accuracy independently of particle content. 



Henceforth energy density must be understood as that in 
the cluster's rest frame, that is the ratio between mass and 
proper volume M/V . The standard value of 0.4 GeV/fm 3 has 
been chosen and used throughout this paper as it corresponds, 
in the thermodynamical limit, to the energy density of a hadron 
gas at a temperature of about 160 MeV, which has been de- 
termined in previous analyses of particle production in high 
energy collisions |11) . 
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The actual CPU time 3 needed for the calculation of 
f2{Nj} with 1000 random samples is reasonably small. This 
is shown in fig. [21 for channels with N neutral pions in a 
cluster with mass JV(m T o + 500 MeV) and energy density 
0.4 GeV/fm 3 . In a full quantum statistics calculation, in- 
cluding lumps up to 5 particles, the CPU time is much 
larger because of the extra terms in the cluster decompo- 
sition. The time needed to calculate the integral J7J in- 
creases almost exponentially with the number of particles 
N at low N, whereas at N = 11 it features a drop due to 
a switch in the numerical method (see Appendix A) and 
flattens out thereafter. In quantum statistics, even at high 
N, there are many phase space integrals in the cluster de- 
composition to be calculated with the same method as for 
the classical terms at low N, thus the CPU time steadily 
increases, though not as steeply as for the classical term 
alone. 



4 Importance sampling of microcanonical 
ensemble 

As has been discussed previously, our goal is to develop a 
fast and accurate numerical method to calculate averages 
like l|10|) in the microcanonical ensemble. Being able to ef- 
fectively calculate ^{n,} for any channel, a brute force op- 
tion is to do it for all of them. However, this method is not 
appropriate for a system like the hadron gas, because the 
actual number of channels is huge. Indeed, with 271 light- 
flavoured hadrons and resonances (those included in the 
latest Particle Data Book issue fTJ]), the number of chan- 
nels allowed by energy-momentum conservation is enor- 
mous and it increases almost exponentially with cluster 
mass (see fig. 01, involving an unacceptably large comput- 
ing time. For instance, the CPU time needed to compute 
^{Nj } f° r all of the 23 millions of channels of a cluster with 
4 GeV mass is around 400 hours. Charge constraints can 
indeed reduce significantly the number of allowed chan- 
nels, yet not enough. Therefore, the calculation of the 
phase space volume of all allowed channels is possible only 
for very light clusters, in practice lighter than ss 2 GeV. 
Hence, if a method based on the exhaustive exploration 
of the channel space is not affordable, one has to resort 
to Monte-Carlo methods, whereby the channel space is 
randomly sampled. 

An estimate of the average (|10|) can be made by means 
of the so-called importance sampling method. The idea of 
this method is to sample the channel space (i.e. the set of 
integers Nj, one for each hadron species) not uniformely, 
rather according to an auxiliary distribution iZ^jy.} which 
must be suitable to being sampled very efficiently to keep 
computing time low, and, at the same time, as similar as 
possible to the distribution f2^ Nj y. The latter requirement 
is dictated by the fact that /2{jv-} is sizeable over a very 
small portion of the whole channel space. Hence, if ran- 
dom configurations were generated uniformely, for almost 

3 All CPU times quoted throughout are referred to a personal 
computer with Pentium IV processor working at 2 GHz clock 
rate, with rated SPECint=426 and SPECfp=304 




M (GeV) 

Fig. 3. Number of allowed channels as a function of cluster 
mass in a hadron gas with 271 species and free charges. 



all of them f2{jsr y would have a negligible value, thus a 
huge number of samples would be required to achieve a 
good accuracy. On the other hand, if samples are drawn 
according to a distribution similar to J?{ jy } ; little time is 
wasted to explore unimportant regions and the estimation 
of the average (|1()|> is more accurate. A crucial requirement 
for is not to be vanishing or far smaller than Q{Nj} 

anywhere in its domain in order not to exclude some good 
regions from being sampled, thereby biasing the calculated 
averages in a finite statistics calculation. 
Rewriting Eq. (fT77|) as: 



(O) 



n 



{Nj} 



E 



{Nj} 

• ; v ; 77^77 //;v QN; v<1 



(16) 



makes it apparent that a Monte-Carlo estimate of (O) is: 



i(*0 



Q 



(fe) 



(O) 



(fe) 

{Nj} 



(17) 



n 



(fe) 

{Nj} 



where {Nj}^ are samples of the channel space extracted 
according to the distribution 77 and fulfilling the charge 
constraint Q = J^j Njqj. 

Provided that M is large enough so that the distribu- 
tions of both numerator and denominator in Eq. I|17f) are 
gaussians (hence the conditions of validity of the central 
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limit thoerem are met), the statistical error <J(o) on the 
average (O) can be estimated as (see Appendix B): 



1 



(Of 



-n 



O 



2 n {N 3 } 
n {N 3 }, 



(0) 2 J? 2 



-2(0) 



-n 



-n 



Q {N 3 } 



o 



n {Nj}. 



- (0}f2 2 



(18) 



where E77 stands for the expectation value relevant to the 



77 distribution. If Q 



n 



reduces to the familiar form 
1 



the above expression 



a 



(O) 



M 



(o 2 ) (oy 



(19) 
and 



The estimator (|17|1 has a bias unless Q{n } = -^{jv } 
this confirms the necessity to find a distribution as similar 
as possible to fi^.y. However, the bias, whose general 
expression is derived in Appendix B, scales with 1/M and, 
therefore, for large M, can be neglected with respect to 
the statistical error scaling with \j\[M. 

A possible option for 27/ jv} is the product of K (as 
many as particle species) Poisson distributions: 



A" 



n {N 3 } = n ex p[-^]^T 



(20) 



which will be henceforth referred to as the multi-Poisson 
distribution or MPD. This distribution can indeed be sam- 
pled very efficiently and it is the actual multi-species mul- 
tiplicity distribution in the grand-canonical ensemble, in 
the limit of Boltzmann statistics. Although this distribu- 
tion is not the limit of f2{N } m the thermodynamic limit, 
i.e. for large mass and volume of the clusters ^2] (see 
next section), the similarity between MPD and the corre- 
sponding microcanonical distribution should be sufficient 
to effectively remove the region of the channel space where 
&{Nj} is practically vanishing. The mean values of the 
MPD are free parameters to be set in order to maximize 
the similarity between MPD and 17{jy .}. The most sensi- 
ble choice is to enforce as mean values the mean hadronic 
multiplicities calculated in the grand-canonical ensemble 
with volume and mean energy equal to the volume and 
mass of the cluster. Indeed, unlike the higher order mo- 
ments, the mean multiplicities, or first order moments of 
the multi-species distribution, in the microcanonical en- 
semble converge to the corresponding values in the grand- 
canonical ensemble in the thermodynamical limit, as it 
will be shown later on. These can be calculated, in the 
Boltzmann statistics by a well known formula: 



(2Jj + 1)V 
2tt 2 



(?) n 



A'/ 



(21) 



where V is the cluster's volume, T is the temperature 
and \i the fugacity corresponding to the charge Qi. Tem- 
perature and fugacities are determined by enforcing the 



grand-canonical mean energy and charges to be equal to 
the actual mass M and charges Q of the cluster: 



d 



with 



Zj (T) = 



{2.], + l)V , 



2tt 2 



(22) 



(23) 



The l|22f) are just the saddle-point equations for the asymp- 
totic expansion of the microcanonical partition function, 
showing that the microcanonical ensemble can be approx- 
imated by the grand-canonical ensemble for large masses 
and volumes 

It should be stressed that the l(2"T|) is just a particu- 
lar choice of the za,'s in Eq. I)2()[l. which is by no means 
compelling. The only purpose and merit of the i/j's is to 
make the multiplicity distribution 77 {Nj} m Eq. 1(20(1 as 
close as possible to i^tjy-} to speed up the computation. 
If a choice different from Eq. l]21[l can do a better job, 
this could and should be retained. For the same reason, 
it makes little sense to use the precise canonical expres- 
sions of the z/j's corrected for quantum statistics, because 
at the energy density we are interested in, this is just a 
correction. Altogether, the means in Eq. 1)21(1 turn out to 
be satisfactory for most practical purposes. 

If cluster charges are unconstrained, only the first of 
the equations (122[l is needed and fugacities can be dropped 
(i.e. they are taken to be 1) from Eq. 1)21(1. On the other 
hand, if cluster charges are fixed, among all configurations 
drawn from MPD, only those fulfilling charge conservation 
should be retained and considered for microcanonical av- 
erage calculations in Eq. 1(17(1 . This preselection of random 
samples significantly affects the overall efficiency of the al- 
gorithm, because the acceptance rate of samples extracted 
from MPD becomes small, and it is increasingly smaller 
for larger clusters. The acceptance rate can be improved 
by resorting to a conditional probability decomposition 
technique, described in the next section. 

Summarizing, the procedure to estimate (O) in the im- 
portance sampling method for a given cluster is as follows: 

1. calculate T and Vj's according to Eqs. 1(22 (1 .1(23 (1 and 

2. sample the MPD ((20|) some large number of times and 
for each sample {Nj} compute numerically the inte- 
gral Q{n } by the method described in Sect. 3 with a 
suitable number of Monte-Carlo extractions N$; 

3. calculate the sums in Eq. l(T7|l and estimate the statis- 
tical error according to Eq. 1(18(1 . 

An improvement in the accuracy of the estimation of (O) 
can be obtained by drawing random samples in an ex- 
tended space instead of calculating f2{N } for each chan- 
nel at each step. The idea is to perform the Monte-Carlo 
importance sampling in the variables: 



{JVi 



,N K \r u ...,r N -i} with N = J2 N i ( 24 ) 
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at the same time, instead of calculating f2{N-y separately 
for each extracted sample in {Nj}. We first note that Eq. 
Ijl5(l can be further rewritten as, by using Eq. (|13|l : 



Q 



{Nj} 



E 



n 



( T l) N J+ H i(2Jj + 



i 1I, V I", K i 

yHj,3H-i 



Ah„ 



{h ni },...,{hn K } 

I dri . . . / dr H -i s„ H 

Jo Jo i^) m 

/ dri ... / dr w _! V] 



■T(ri,...,r H _i) 



{hnj,...,^^} 



n 



( T X) JV > +fl >(2J J - + 1)^ 



d n ... f dr^-i ^{N^lin}) 



o 



T(n, . . .,ru-i) 



(25) 



being {r;} = (ri, . . . , rjv-i)- In practice, each term in the 
sum has been multiplied by 1 = dr# ■ ■ ■ Jq drjv-i, with 
H < N and, doing so, we have been able to take out 
the integration on the r^'s. The Monte-Carlo estimate of 
J7/jv } as expressed by the last integral in Eq. (|25|) can 
then be written as: 



{Nj} 



, N S 

— Y*(r<i 



(ft) 



r {k) ) 

' N-l) 



(26) 



fc=i 



(k) 

where r\ are random numbers uniformely distributed in 
the interval [0, 1]. Likewise, looking at Eqs. fTTHl and (J22JI, 
the estimator of the mean value of the observable O in 
this extended sampling space can be rewritten as: 



(O) 



(ft 

{JV 3 -} 



v m' mNj} (k) \{n} {k) ) 

l^k=l 



(27) 



77 



(ft) 

{Nj} 



For a fixed number of calls to the random number gen- 
erator, this method optimizes the accuracy because most 
of them are spent to calculate ^{n,} f° r most probable 
channels, i.e. those for which an improvement in accuracy 
is more rewarding, whereas, in the previous method, ^{jv 3 -} 
was calculated with a hxed number of random samples N$ 
regardless of its size. More specifically, if in the previous 
approach Ns x M computation of the function ^ in Eq. 
H25|) were performed, about the same CPU time is needed 
to calculate M' — Ns x M random samples in Eq. (|27(l . 
but with a sizeable reduction of the statistical error on 
(O). 

In our calculation, resonances whose width exceeds 1 
MeV are handled as particles with a distributed mass. 
For each random sample, the function 'P in Eq. (|25|l is 
calculated by randomly drawing masses from a relativistic 
Brcit-Wigner distribution for each resonance: 



BW(m 2 )dm 2 



1 



(m 2 



ml) 2 



rh 



■ dm 



(28) 



The mass range is symmetric around the central value mo 
with half-width min{2_T, m cu t} where m cut is the minimal 
mass allowed by the known decay channels of the reso- 
nance. If the sum of the extracted masses of particles and 
resonances exceeds the cluster mass, the function ^ is set 
to zero. 



5 Comparison between microcanonical and 
canonical ensemble 

The importance sampling method with MPD provides a 
suitable technique to calculate averages in the microcanon- 
ical ensemble. In this section we will mainly show numer- 
ical results on the difference between averages in the mi- 
crocanonical and canonical ensemble. 

We have calculated average multiplicities in a hadron 
gas including 271 light-flavoured hadron species up to a 
mass of about 1.8 GeV quoted in the 2002 Particle Data 
Book issue f° r completely neutral clusters (Q = 0) 
and pp-like clusters, i.e. with net electric charge Q = 2, net 
baryon number B = 2 and vanishing net strangeness. The 
energy density in the rest frame of the cluster M/V has 
been set to 0.4 GeV/fm 3 , corresponding, in the thermo- 
dynamical limit at vanishing chemical potentials, to the 
temperature value of about 160 MeV found in analyses of 
particle multiplicities in high energy collisions . No ex- 
tra strangeness suppression factor has been used here, as 
we are just interested in a comparison between calculated 
multiplicities in two different ensembles. 

Whilst the energy density was kept fixed, the mass has 
been varied from 2 to 14 GeV for neutral cluster and from 
4 to 14 GeV for pp-like clusters in steps of 2 MeV. For each 
mass, 10 7 random samples fulfilling charge conservation 
(i.e. passing charge preselection) have been drawn from 
the MPD. 

In order to improve the performance of the algorithm 
and decrease the rejection rate at the preselection stage, 
we have implemented a multi-step extraction procedure 
taking advantage of well known properties of the Pois- 
son distribution. Instead of extracting all particle multi- 
plicities independently from Poisson distributions, we ex- 
tracted first the number of baryons Nb and antibaryons 
Ng from two Poisson distributions, with means equal to 
the sums of all baryon and antibaryon means respectively, 
and accept the event only if Nb~N b — B. Indeed, denot- 
ing by TTj a single-species Poisson distribution, the MPD 
constrained with charge conservation can be written as: 

K 



bar antibar mesons 

= ir B (N B )w B {N B )P(Ni,N2, ■ ■ ■ \N B )P(N l ,N- 2 , . . . \N B ) 

X Jj7T(JVj) SNB-N^BSZjNjS^S^NjQja (29) 
mesons 

Here the probability distribution Ilbar ""(-^Y?) °f navm g 
given baryon multiplicities has been decomposed into the 
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product of a Poisson distribution ttb for the overall baryon 
multiplicity Nb and the conditional probability P(Ni, AT 2 , . 
of having the same set of baryon multiplicities provided 
that their sum is Nb', similarly for antibaryons. The P 
distribution is actually a multinomial distribution: 



P(N 1 ,N a ,...\N B )<xN B l*[[^ n 



(30) 



where Vj are given by Eq. (|21|) . Once all baryon and an- 
tibaryon multiplicities have been extracted according to 
the multinomial distributions, which can be sampled as 
quickly as Poisson's, the same procedure is carried out for 
strange mesons. Since strange and antistrange mesons are 
independent of the previous baryon extraction, the num- 
ber of strange mesons Ng and antistrange mesons N§ are 
extracted according to a Poisson distribution subject to 
the requirement that N$ — Ng — S — Sb, where Sb is 
the net overall strangeness carried by the previously ex- 
tracted baryons. If Ns and N§ fulfill the above condition, 
strange and antistrange meson multiplicities are extracted 
from a multinomial distribution, otherwise the event is re- 
jected and the whole extraction procedure gets back to the 
very beginning, i.e. by randomly sampling baryon and an- 
tibaryon numbers. Likewise, the number of charged non- 
strange mesons and their antiparticles is extracted and, 
provided that charge conservation is fulfilled, their indi- 
vidual multiplicities are determined. Finally, the multi- 
plicities of completely neutral mesons, which do not affect 
the total charges, are extracted. The advantage of this 
method over that based on straight multi-Poisson sam- 
pling can be more easily understood in terms of rejection 
rate. In fact, the number of rejected channels {Nj} out of 
those sampled because of the charge constraints, should 
be the same in both methods as they are drawn from the 
same distribution, i.e. the MPD. Nevertheless, in the lat- 
ter multi-step method, a single rejection of a channel on 
average does not require K (as many as hadron species) 
extraction from a Poisson distribution: it may well occur 
at the first stage of baryon number conservation, with only 
one extraction from a Poisson distribution or later, with 
a number of sampled Poisson or multinomial distribution 
which is still less than K. The actual gain in CPU time 
may be dramatic, especially for high mass clusters, where 
the number of allowed channels is huge. For instance, for 
a pp-like cluster of 14 GeV, we have estimated a ratio of 
0.025 between the average CPU time needed to extract a 
channel with proper charges in the multi-step improved 
method and in the original method. 

Altogether, the actual CPU time needed to generate 
10 7 accepted samples (i.e. the statistics relevant to the 
plots in figs. 1415( 1 of multi- hadronic channels and calculate 
average multiplicities in the importance sampling method 
for a neutral 4 GeV cluster at 0.4 GeV/fm 3 energy density 
is about 4.6 10 2 sec, i.e. less than 8 minutes. 

The canonical average multiplicities to be compared 
with the microcanonical ones have been calculated by de- 
termining first the temperature corresponding to the sad- 
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< 
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Fig. 4. Relative difference between microcanonical and canon- 
ical average primary multiplicities of mesons ((-Wj) micro — 
{Nj) C an)/{Nj) C an for completely neutral clusters with mass 4, 
8, 12 GeV . The error bars indicate the statistical error of 
the importance sampling Monte-Carlo computation. Connect- 
ing lines are drawn to guide to eye. 



die point equation 0: 



d 



M-T 2 — log Z(Q,T) = 



(31) 



+ 7T 



d 3 </> e 1 ^ exp 



(2Jj + l)V 



x /d 3 p log(l±e-VP a+m 3 -/ T -"fc-*) ±1 l (32) 



where 
Z(Q,T) 



is the canonical partition function; as usual, the upper 
sign applies to fermions, the lower to bosons. Then, mul- 
tiplicities have been calculated according to the known 
expression in the canonical ensemble |3]: 

(N,) = 

n— 1 v y 

where terms in the series beyond n = 1 account for quan- 
tum statistics effects and are in fact important only for 
pions at the usually found temperature values of 160- 
180 MeV. As for the microcanonical ensemble, resonances 
whose width exceeds 1 MeV have been considered as free 
hadrons with a mass distributed according to a relativistic 
Breit-Wigner distribution. 
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Fig. 5. Relative difference between microcanonical and canon- 
ical average primary multiplicities of baryons ((-/Vj) m i cro — 
{Nj) ca , n )/ (Nj) ca _ n for completely neutral clusters with mass 4, 
8, 12 GeV. The error bars indicate the statistical error of the 
importance sampling Monte-Carlo computation. Connecting 
lines are drawn to guide to eye. 



Fig. 7. Relative difference between microcanonical and canon- 
ical average primary multiplicities of baryons ((iV,) m i Crc , — 
(Nj) can) I (JVj)can for pp-like clusters with mass 4, 8, 12 GeV. 
The error bars indicate the statistical error of the impor- 
tance sampling Monte-Carlo computation. Connecting lines 
are drawn to guide to eye. 
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Fig. 6. Relative difference between microcanonical and canon- 
ical average primary multiplicities of mesons ((A f j} m i cro — 
(Nj) can) /{Nj) can for pp-like clusters with mass 4, 8, 12 GeV. 
The error bars indicate the statistical error of the impor- 
tance sampling Monte-Carlo computation. Connecting lines 
are drawn to guide to eye. 



Fig. 8. Relative difference between microcanonical and canon- 
ical average primary multiplicities of antibaryons ((iVj) m i Cro — 
{Nj)ca.n)/{Nj)ca.n for pp-like clusters with mass 4, 8, 12 GeV. 
The error bars indicate the statistical error of the impor- 
tance sampling Monte-Carlo computation. Connecting lines 
are drawn to guide to eye. 
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Fig. 9. Relative difference between microcanonical and canon- 
ical average primary multiplicities of it , p and fl ((iVj) micro — 
{Nj) ca , n )/{Nj) ca _ n for a neutral cluster with mass of 6 GeV as 
a function of energy density. The error bars indicate the sta- 
tistical error of the importance sampling Monte-Carlo compu- 
tation. Connecting lines are drawn to guide to eye. 
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Fig. 10. Relative difference between microcanonical and 
canonical average primary multiplicities of 7r° ((A^) m i cro — 
(Nj} ca , n )/(Nj} can a completely neutral cluster and in a pion 
function of cluster mass at an energy density of 0.4 
GeV/fm 3 . The error bars indicate the statistical error of the 
importance sampling Monte-Carlo computation. Connecting 
lines are drawn to guide to eye. 



The relative difference between microcanonical and ca- 
nonical average for hadrons belonging to basic SU(3) mul- 
tiplcts and three different cluster masses are shown in figs. 
I4I5I6I7I and IH1 At a mass of 4 GeV, the deviation is signif- 
icant and shows a considerable variation as a function of 
the species. In neutral clusters, mesons and baryons fea- 
ture a different behaviour: whilst microcanonical multi- 
plicities of mesons are higher than the corresponding mul- 
tiplicity in the canonical ensemble, those of baryons are 
lower. In pp-like cluster, on the other hand, the general 
behaviour is not as simple; different mesons show differ- 
ent signs of the difference and there are even oscillations 
of it going from low to high mass clusters. On the other 
hand, it is evident that already at 8 GeV of mass, where 
the total primary multiplicity of particles is around 8 (see 
table 2), all differences between the microcanonical and 
canonical ensemble do not generally exceed 20% and fur- 
ther shrink to about 10% at 12 GeV. Furthermore, these 
differences depend only weakly on the energy density, as 
shown in fig.|5| We can conclude that, as a rule of thumb, 
the canonical ensemble is a good approximation of the 
microcanonical one for masses > 8 GeV at energy densi- 
ties between 0.1 and 0.9 GeV/fm 3 with quantum numbers 
not greater than that of an elementary colliding system. 
Provided that the further assumption of the equivalence 
between the actual set of clusters and an equivalent global 
cluster (EGC) holds [B], this result justifies the use of the 
canonical ensemble in the analysis of particle multiplici- 
ties in pp, e + e _ and other high energy collisions. Indeed, 



the mean masses of the EGC corresponding to the actu- 
ally fitted temperatures and volumes, shown in tabic ^ 
turn out to be sufficiently large for most collisions, with 
the likely exception of K+p at y/s = 11.5 GeV and e + e~at 
^fs = 14 GeV where the mean mass is lower than 7 GeV. 
At y/s > 20 GeV, the mean mass is larger than 8 GeV, 
implying that the canonical ensemble is a good approxi- 
mation. 

The quick convergence of microcanonical average mul- 
tiplicities to canonical ones in a hadron gas is favoured 
by the large number of available degrees of freedom. From 
a mathematical point of view, a large number of degrees 
of freedom makes the saddle point expansion converging 
faster. In physical terms, there are more ways to conserve 
energy and momentum in a system with a larger num- 
ber of particle species, so that fulfilling these constraints 
become less important earlier than in a system with few 
degrees of freedom. This is demonstrated in fig. 1101 where 
the relative difference between neutral pion multiplicity 
in microcanonical and canonical ensemble is shown for a 
completely neutral pion gas and for a full hadron-g as as a 
function of cluster mass. 

The microcanonical average multiplicities of most had- 
rons, to a good approximation, increase linearly as a func- 
tion of cluster mass, for fixed energy density, starting from 
M ~ 3 GeV. On the other hand, particles with large 
charge content (such as fl) show a stronger dependence 
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Table 1. Fitted parameters temperature, volume and extra strangeness suppression parameter (7s or the mean value of 
produced strange quark pairs out of the vacuum (ss)) in elementary collisions in the canonical analysis of hadron abundances. 
The canonical fit assumes the equivalence, as far as particle multiplicities are concerned, between the set of actual clusters and 
one global cluster whose resulting mean mass is quoted in the last column. Note that e + e~ — > cc and e + e~ — > bb events, 
where heavy flavoured hadrons are emitted, have been excluded to estimate (M) in e + e~ collisions. 



Collision 


(GeV) 


Reference 


T (MeV) 


V (fm 3 ) 


7S or (ss) 


(M) (GeV) 


K+p 


11.5 




6 




176.9±2.6 


8.12±0.83 


(ss)= 0.347±0.020 


6.51 


K+p 


21.7 




6 




175.8±5.6 


12.0±2.4 


(ss)= 0.578±0.056 


8.47 


7T + p 


21.7 




6 




170.5±5.2 


16.7±3.1 


(ss)= 0.734±0.049 


8.23 


PP 


17.2 




14 




187.2±6.1 


6.79±1.6 


(ss)= 0.381±0.021 


7.74 


PP 


27.4 




6 




162.4±5.6 


25.5±1.8 


(ss)= 0.653±0.017 


9.67 


e + e~ 


14 




6 




167.3±6.5 


15.9±4.1 


7s= 0.795±0.088 


6.08 


e + e~ 


22 




6 




172.5±6.7 


15.9±4.8 


7s= 0.767±0.094 


8.12 


e + e~ 


29 




6 




159.0±2.6 


33.1±4.0 


7s= 0.710±0.047 


9.28 


e + e~ 


35 




6 




158.7±3.4 


33.7±5.2 


7s= 0.746±0.040 


9.54 


e + e~ 


43 




6 




162.5±8.1 


29.0±9.2 


7s= 0.768±0.065 


9.99 


e + e~ 


91.25 




11 




159.4±0.8 


52.4±2.2 


7s= 0.664±0.014 


16.0 


PP 


200 




11 




175±11 


35±14 


7s= 0.491±0.044 


21.6 


PP 


546 




11 




167±11 


65±27 


7s= 0.526±0.044 


28.7 


PP 


900 




11 
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Fig. 11. Microcanonical average primary multiplicities of 7T , 
p, n and of all hadron species (h) in a completely neutral 
cluster as a function of cluster mass at an energy density of 
0.4 GeV/fm 3 . The error bars indicate the statistical error of 
the importance sampling Monte-Carlo computation. Lines are 
drawn to guide to eye. 



on mass, as shown in fig. ^] over the mass range appro- 
priate for microcanonical and canonical calculations. 

We have also compared the overall primary multiplic- 
ity distributions in the two ensembles. The multiplicity 
distribution in the microcanonical ensemble has been de- 



termined by taking Na as observable in Eqs. tfTT7f> . (fTfi|) 

and performing an importance sampling Monte-Carlo cal- 
culation. The multiplicity distribution in the canonical en- 
semble has been determined by the same method, that is 
extracting particle numbers from the MPD and weigh- 
ing each event fulfilling charge conservation with the ratio 
between the actual multi-species multiplicity distribution 
and the MPD. The former reads (see Appendix C): 



P({Nj}) 



1 



Z(Q) 



n e ^ Nj+Hi 
3 {h nj } 



n 



=1 n. 



'h I 



v., (34) 



where {h nj }, as usual, denote partitions and: 



z i{n s ) 



(2tt) 3 



P e 



(35) 



For T w 160 MeV, only the leading poissonian terms in 
the distribution l|34() corresponding to {h nj } — (N,0, . . .) 
can be retained for all particles except pions. 

The comparison between the multiplicity distributions 
are shown in figs . 1 1 21 and 1 1 31 for neutral and pp-like clus- 
ters respectively. It can be seen that the mean values tend 
to get closer as mass increases, whilst the dispersion is 
lower in the microcanonical than in the canonical ensem- 
ble, where the distribution is almost poissonian. This re- 
markable effect on particle number fluctuation is owing to 
the overall energy-momentum constraint causing a global 
correlation in particle production. The ratio of the disper- 
sion to the square root of the mean (i.e. the dispersion 
of a Poisson distribution) tends to a factor 1/2 in the 
thermodynamical limit (see fig. I14J) . thus showing that 
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Fig. 12. Comparison between microcanonical (solid) and 
canonical (dashed) overall primary multiplicities distribution 
in neutral clusters of four different masses at an energy den- 
sity of 0.4 GeV/fm 3 . 



Fig. 14. Ratio between the dispersion of the microcanonical 
multiplicity distribution and the square root of its mean in 
neutral (top) and pp-like (bottom) clusters as a function of 
mass at an energy density of 0.4 GeV/fm 3 . Lines are drawn to 
guide the eye. 
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Fig. 13. Comparison between microcanonical (solid) and 
canonical (dashed) overall primary multiplicities distribution 
in a pp-like cluster of four different masses at an energy den- 
sity of 0.4 GeV/fm 3 . 



microcanonical ensemble is not equivalent to the grand- 
canonical ensemble with respect to particle number fluc- 
tuations. We are not currently aware of a simple reason of 
this fact. 

Because of the persistence of the shape difference be- 
tween multiplicity distributions in the two ensembles, the 
number of channels sampled with the MPD whose weight 
is much larger than their corresponding microcanonical 
one increases with total mean multiplicity, hence with 
cluster's mass (at fixed energy density). This is the rea- 
son of the slightly increasing statistical error on average 
multiplicities and multiplicity distributions for increasing 
cluster mass seen in figs. I4I5I6I7I8I12I13I One could rem- 
edy this and speed up microcanonical calculations for very 
large clusters by changing the sampling distribution. For 
instance, instead of using a MPD to sample individual 
species numbers, the total number of particles can be first 
drawn from a gaussian distribution with variance 1/4 of 
the sum of all hadron primary multiplicities as estimated 
in the canonical ensemble; then, the number of baryons 
Nb, strange mesons Ns, charged non-strange mesons Nq, 
neutral mesons TYo and their respective number of antipar- 
ticles could be sampled from a multinomial distribution: 

p = m Uh (36) 

where i = B, B, S, S, Q, Q, and are the sums of all Vj 
functions in Eq. (|21|) of particles belonging to the same 
class i. Finally, once the multiplicities of each class TV* 
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fulfilling conservation laws have been extracted, individual 
particle multiplicities could be determined for each class 
in turn by using again multinomial distributions. 

Table 2. Mean multiplicity (N) and dispersion D of the mi- 
crocanonical and canonical multiplicity distribution for neutral 
and pp-like clusters as a function of the mass. Also quoted the 
corresponding temperature in the canonical ensemble obtained 
from Eq. J3TJ. 



Neutral cluster 





Microcanonical 


Canonical 


M (GeV) 


(N) 


D 


(TV) 


D 


T (MeV) 


2 


2.43 


0.63 


1.96 


1.74 


175.3 


4 


4.45 


1.01 


4.08 


2.45 


169.3 


6 


6.53 


1.24 


6.17 


2.93 


166.3 


8 


8.61 


1.43 


8.27 


3.29 


164.6 


10 


10.67 


1.59 


10.38 


3.58 


163.5 


12 


12.76 


1.75 


12.46 


3.90 


162.8 


14 


14.84 


1.89 


14.56 


4.15 


162.2 


pp-like cluster 




Microcanonical 


Canonical 


M (GeV) 


(TV) 


D 


<7V> 


D 


T (MeV) 


4 


3.87 


0.81 


3.63 


1.33 


142.1 


6 


6.04 


1.08 


5.76 


2.05 


152.5 


8 


8.20 


1.30 


7.90 


2.59 


156.1 


10 


10.35 


1.49 


10.04 


3.02 


157.8 


12 


12.48 


1.66 


12.17 


3.40 


158.7 


14 


14.60 


1.82 


14.27 


3.73 


159.2 



6 The Metropolis algorithm 

The importance sampling method allows to calculate av- 
erages like (|10H in single Monte-Carlo runs quite straight- 
forwardly and can be used as an event generator if events 
are reweighted, as we have seen, by a factor fi^^ / Il^.y 
If one needs to sample J?{ ^ j j. directly without reweighting 
the events, different methods should be considered. A hrst 
possibility is a rejection method, but it can be soon real- 
ized that it is unfit for the problem we are dealing with. 
In fact, for this method to be effective, one needs a cover- 
ing function F({Nj}) (i.e. a function of the channel such 
that F({Nj}) > ft{Nj\ V{A^} which can be sampled very 
efficiently and as close as possible to the J?{jy.}. Such a 
function is hard to find because ^{Nj} is not smooth in its 
domain; strong variations of the phase space volume may 
occur by changing just one particle. We have seen that the 
MPD is likely to be similar to our target distribution, but 
in order to be a covering function, it should be rescaled by 
a factor c such that cF[{ N .} > ^{Nj} V{iVj}. However, to 
estimate c, one ought to evaluate ^{n,} f° r a U channels 
and this is just what cannot be afforded. A general method 
to sample complex multi-dimensional distributions is the 
Metropolis algorithm [TJ], which has been applied to the 
specific problem of numerical calculations of the multi- 



hadronic microcanonical ensemble by Werner and Aiche- 
lin [7] . In this work, we will show that this method can be 
further improved in speed and accuracy and, at the same 
time, that great care is needed in handling it, especially 
when assessing the equilibrium conditions. 

The Metropolis algorithm prescribes the implementa- 
tion of a random walk in the channel space on the basis of 
acceptance or rejection of proposed transitions from the 
current position. After some number of steps, the proba- 
bility of visiting a given point is proportional to the tar- 
get distribution; otherwise stated, the points in the ran- 
dom walk are actual samples of the the target distribution 
(in our case fi{ N .}) and they can be stored as generated 
events. We will now discuss more in detail how this comes 
about. In a general random walk, the probability P m (t) of 
visiting a state m (i.e. a channel or multi-hadronic con- 
figuration) at the t th step evolves according to the master 
equation: 

P m (t+1)-P m (t) = Pn(t)w{n -» m)-P m (t)w(m -» n) 

n 

(37) 

where w(n — > to) is the probability of transition from the 
state n to the state to. At the equilibrium: 

P m (t + 1) = P m (t) Vto (38) 

This occurs if (a sufficient but not necessary condition): 

P n (t)w(n — ► m) = P m (t)w(m — > n) Wi, m (39) 

If we want the probabilities P m to be proportional to 
a given target distribution f(m), so that the number of 
times the state to is visited from some step t onwards is 
proportional to f(m), we have to enforce the condition: 

w^rn) (4Q) 
uu[m —> n) J{n) 

Provided that the above equation is fulfilled, the choice 
of a set of transition probabilities w (a so-called updating 
rule), is free. The only requirement is that every point can 
be reached from any point in a finite number of steps with 
non- vanishing probability (ergodicity condition) , otherwise 
there would be unaccessible regions even though / ^ 
therein. Although any choice of the w's is in principle al- 
lowed, some are worthier. Indeed, the transition probabili- 
ties w govern the dynamical behaviour of the random walk 
and, particularly, how fast the system gets to the equilib- 
rium condition (|38J) after a transient. This is owing to the 
fact that, in general, the random walk starts from states 
which are not random samples of the target distribution, 
so that P m (0) 7^ f{m). A good choice of the w's will keep 
the relaxation time T re i, defined as the number of steps 
needed to get sufficiently close to the equilibrium value, 
to a minimum, thereby making event generation faster. 

In general, the transition probability can be decom- 
posed into a proposal probability T (proposal matrix), i.e. 
the probability of considering a given transition, and the 
conditional probability A of accepting it once it has been 
proposed. In symbols: 

w(n — > to) = T(n — > m)A(n — > to) . 
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Starting from a non-equilibrium situation, a physical sys- 
tem reaches equilibrium earlier if transition probabilities 
are larger. Likewise, in the Metropolis algorithm, the re- 
laxation time is little if w(n — » to) are large; in other 
words, if w(n — > n) is small taking into account the nor- 
malization condition: 



w(n — > to) = 1 



(41) 



For the transition probabilities to states different from the 
current state to be large, the acceptance matrix should be 
as large as possible for m ^ n once the proposal matrix 
is known. There is indeed an optimal choice for A which 
reads: 



A(n — > to) = min < 1 



f(m)T(m — > n) 
/(n)T(n — > to) 



(42) 



This choice maximizes w(n — > to) if n ^ m |16| and ful- 
fills the condition Q4Up. The next problem is choose a good 
proposal matrix. The issue is discussed in detail in ref. 
and the conclusion is that T(n — > m) ~ /(to), i.e. the pro- 
posal matrix should be an easy-to-sample distribution as 
close as possible to the target distribution /(to). This can 
be easily understood by taking the limiting case (which, 
if possible, would make the Metropolis algorithm unnec- 
essary) : 

T(n ->m) = /(to) (43) 

In this case, the relaxation time would be zero: the start- 
ing point, as well as all other points in the random walk, 
would be sampled from the target distribution itself and 
any transition would be accepted because A(n — > to) = 
1 according to Eqs. and (@2J>- While the condition 
cannot be obtained in practice (Metropolis algorithm 
would be unnecessary in that case), one can try to get 
as close as possible to it. From what we have seen in 
Sect. 4 about importance sampling, it can be argued that 
the MPD in Eq. 1)20(1 would be a good proposal matrix 
and could also be used to pick the starting point in the 
Metropolis random walk. We will see the benefits of the 
use of the MPD more in detail in the next section. 

A drawback of the Metropolis algorithm is that, un- 
like in the importance sampling method, different samples 
(i.e. steps in the random walk) are not independent. This 
can be easily understood reminding that proposed transi- 
tions may be rejected, so that a point may appear several 
times in a row in the random walk. Hence, there is a fi- 
nite positive statistical correlation between the values of 
physical observables at different steps and this gives rise 
to an increase of the overall uncertainty in the estimate of 
averages with respect to the case of independent samples. 
The fact that different events are correlated may render 
the use of Metropolis algorithm not appropriate in some 
applications. However, in the problem of high energy col- 
lision event simulation, where one has to hadronize many 
clusters with different masses in each event, this is not an 
issue; in this case, a single Metropolis random walk must 
be run for each cluster, and only one sample, representa- 
tive of its microcanonical ensemble, drawn after equilib- 
rium has been reached. 



The Metropolis algorithm can be used to estimate the 
average of an observable like l|l()(l in the microcanonical 
ensemble of a single cluster by taking a sufficiently large 
number of steps (i.e. 3> T re {) in one random walk and 
calculating: 



(O) 



M 



(44) 



where M is the total number of steps and 0^ k ' the actual 
value of the observable O at k th step. A nice feature of 
the Metropolis algorithm is that, unlike in the importance 
sampling method, there is no need of overall normalization 
when estimating (O) (compare Eq. (|44|) with Eq. I|17|)l. On 
the other hand, as already emphasized, the values of O at 
different steps are correlated. The statistical error on (O) 
in Eq. I|44J) can be estimated as (see Appendix D): 



r (o> 



<o 2 ) - (oy 



2R 



(45) 



where R is the integral of the autocorrelation function 
A, defined here as the difference between the expectation 
value of the product of the observable value at the steps 
k and k + h and its average (O) 2 



A(h) = E(0 (,£) (fc+ ' l) ) - {Of 



(46) 



If the starting point is random, the autocorrelation func- 
tion is independent of k and gives information about how 
correlated are distant steps, namely how long the system 
keeps memory of its past steps. As already mentioned, 
this correlation between different steps arises from the fi- 
nite probability of rejecting transitions and vanishes only 
if all proposed transitions are accepted, i.e. if Eq. (|43|l is 
fulfilled. Therefore, the autocorrelation function is always 
positive and vanishes only if different steps are indepen- 
dent. Moreover, the statistical error (|45|l is larger than 
in the case of uncorrelated samples, where it reaches its 
minimum. The number of steps needed to reduce A(h) to 
some small fraction of (O) 2 is defined as autocorrelation 
time T auto - Thus, in order to minimize the statistical error 
on (O) in Eq. I|45l) . one should keep the autocorrelation 
time as little as possible. 

The autocorrelation function and time depend only on 
the updating rule, whilst the relaxation time T re i also de- 
pends on how the initial state is chosen. In principle, the 
relaxation time might be zero while the autocorrelation 
time is not. When using the MPD both as proposal matrix 
and to generate the starting point, these two quantities are 
tightly related. 

The autocorrelation function can be estimated during 
the Metropolis random walk by the sum: 



A(h) 



Er=i uto o^o^+v i 



M 



M - T„ 



\fc=i / 



(47) 



and it is shown in fig. ^] for a cluster of 4 GeV with 
updating rule based on the MPD. 

Now the question arises whether the Metropolis al- 
gorithm leads to more or less accurate computations of 
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Fig. 15. Autocorrelation function for the total primary multi- 
plicity for a neutral cluster of 4 GeV mass at an energy density 
of 0.4 GeV/fm 3 . 

mean values than the importance sampling method, for 
given computing resources and using the same distribu- 
tion Il^jsr y as proposal matrix and sampling distribution 
respectively. This is studied in more detail in Subsect. 7.4. 



7 Study of the Metropolis algorithm 

We have studied the capability of the Metropolis algo- 
rithm as an event generator for the statistical model of 
hadronization and as a computing tool for the hadron gas 
microcanonical ensemble. As has been discussed in the 
previous section, the benchmark for this algorithm is the 
number of steps needed to reach equilibrium, or relaxation 
time T re i , which must be kept to a minimum so as to draw 
samples of the target distribution as quickly as possible. 
Likewise, the autocorrelation time T auto must be small in 
order to minimize the statistical error Ij45(l on the estimate 
of mean values in single random walks. As has been men- 
tioned at the end of previous section, these two quantities 
are tightly related in our MPD-based scheme, so we can 
confine ourselves to study T re i . 

The relaxation time in principle depends on the phys- 
ical observable (average multiplicities, multiplicity distri- 
bution etc.) and it is not easy to estimate in advance. 
Hence, in practice, it must be determined a posteriori 
by analyzing the convergence to equilibrium for the ob- 
servable of interest. For instance, if we ought to calculate 
the overall multiplicity distribution P n , we would have to 
study the height of the n th bin, for all n's, as a function 
of the step. The height of the n th bin at the k th step, that 
is P n (k)i can be estimated by repeating the Metropolis 
random walk many times and averaging: 

! L 

Pn{k) = ^Z)*n,n iW (48) 
i=l 

where rij(M is the actual multiplicity in the i th random 
walk at the step k. 
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Fig. 16. Step hystogram showing the convergence to equilib- 
rium of the probability P7 of a channel with 7 primary hadrons 
as a function of the step in a set of 500000 Metropolis random 
walks for a neutral cluster of 6 GeV with energy density of 0.4 
GeV/fm 3 . The horizontal solid line indicates the P7 value esti- 
mated with the importance sampling method and the dashed 
band its relevant statistical uncertainty. The arrow points to 
the equilibrium point determined through the WOSSR test 
(see text). 



Establishing when a stable value of the examined ob- 
servable is attained can be done by studying its average 
over many Metropolis random walks, like l|48|l as a func- 
tion of the step, i.e. forming a hystogram (the step hys- 
togram, see fig. I16H with Metropolis step numbers as bins 
and the average of the observable (e.g. multiplicity) as 
bin content. In all considered cases, the step hystogram 
shows the same general behaviour, namely a tendency to 
an equilibrium value after, possibly, a strong initial os- 
cillation and followed, sometimes, by very mild damped 
oscillations. Step hystograms show fluctuations around an 
apparent equilibrium value at some scale, as it can be 
seen in figs. 1161171 The fluctuation pattern is determined 
by the interplay of finite statistics and possible dynami- 
cal oscillations governed by the master equation l)37|l . The 
amplitude of statistical fluctuations is a function of the 
updating rule. 

Because of these facts, providing a good quantitative 
assessment of where stability is achieved is not straight- 
forward. Were the statistics of Metropolis random walks 
infinite, the criterium for stability would be fully deter- 
ministic, that is based on the relative difference between 
asymptotic value (the true microcanonical average) and 
actual value from a given step onwards. On the other hand, 
we can only deal with finite statistics, so the estimation 
of a stability point (thence T re i) for a given observable 
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can be done statistically, thus it will be affected by some 
uncertainty. 

The statistical method we have chosen in order to as- 
sess the convergence to equilibrium is based on a non- 
parametric statistical test, the Wilcoxon one sample signed 
rank test (WOSSR), that we will shortly describe in the 
following. The WOSSR test is a non-parametric test, i.e. 
it does not require the knowledge of the statistical distri- 
bution of the sample, and the hypothesis under test (the 
null hypothesis) is that a set of M independent random 
samples {x%, . . . ,xm} has a given median X m . The only 
requirement is that their distribution is symmetric. The 
test procedure consists in sorting the absolute values of 
the differences di = x% — X m \/i and giving them a rank, 
namely 1 to the largest, 2 to the second largest and so 
on. If X m is a good median, the sum of the ranks of the 
positive differences is expected to be close to the sum of 
the ranks of negative differences. The test statistic is just 
the sum W of the ranks of the positive differences and the 
test result is good when W is close to half the sum of the 
first M integers, that is M(M + l)/4. 

We have used this method to study the convergence to 
equilibrium in the Metropolis algorithm taking as random 
samples a subset of M running bins (from the k th to the 
(k + M ) th ) in the step hystogram of the multiplicity distri- 
bution P n (k) (see fig. 116(1 . In fact, if equilibrium is reached 
in the Metropolis random walk (as it is apparent in the 
rightmost part of the hystogram in fig. Ilfijl . P n (k) is ex- 
pected to evenly fluctuate around the asymptotic equilib- 
rium value, whereas a net drift towards this value appears 
when out of equilibrium (i.e. in the leftmost part of fig.H6|). 
Therefore, provided that the distribution of fluctuations is 
symmetric at equilibrium, the asymptotic value is likely to 
be a good median for sets of P n nA with k = 1, . . . , M in 
the equilibrium region and not in the drifting region. Ac- 
cordingly, the WOSSR test will yield a good confidence 
level in the former case and a very small one in the latter. 
The true asymptotic value is not known in practice and 
must be estimated from the step hystogram itself. A good 
choice is the arithmetic mean of a set of hystogram con- 
tents in the rightmost region, where stability is apparently 
reached. 

The non-parametric WOSSR test carried out on a set 
of M running bins in the step hystogram seems suitable 
for our problem. In fact, we do not know the statistical 
distribution of the random variable 0^ (i.e. the exam- 
ined observable) at each step for a given number L of ran- 
dom walks. On the other hand, it is reasonable to assume 
that it is symmetric at equilibrium around the asymp- 
totic value which is one of the requirement of the WOSSR 
test. However, this test also requires the independence of 
samples, which is not the case here because subsequent 
steps in Metropolis random walks are indeed correlated, 
as has been emphasized. Thus, the results of WOSSR test 
in this context should be taken with much care and could 
be misleading if M <C T auto . If a large fluctuation from 
equilibrium value occurs at some point, its persistency in 
sign is fed by the positive correlations of adjacent points 
and this will drive the test to failure even if equilibrium 



was actually achieved. Conversely, it is quite unlikely that 
the test yields a positive answer on a sufficiently large 
set of running bins if equilibrium is not achieved, unless 
two accidentally large fluctuations of equal time size arise. 
Therefore, albeit not appropriate in principle, the WOSSR 
test with M ~ T auto provides a fairly good indication of 
equilibrium. 

In order to actually define an equilibrium point, and 
thereby a relaxation time T re i, we first calculate the arith- 
metic mean of the rightmost 50 (which is of the same order 
as of T auto ) bin contents, which is taken as median X rn to 
be fed in the test. The WOSSR test is then performed on 
subsets of M = 10 running bins in the step hystogram 
taking as starting bin the leftmost and moving rightward 
by one bin at a time. The test is regarded as successful 
if it yields a confidence level of at least 0.05 for a median 
differing from the previously defined X m at most by 0.5%. 
If the test is successful for 10 starting bins in a row, the 
first of those ten is taken to be the equilibrium point. Lit- 
tle variation of the equilibrium point is found by changing 
the number of running bins from 5 to 20. 

Henceforth, we will use the observable overall multi- 
plicity and WOSSR test to point out some important fea- 
tures of the Metropolis algorithm. 

7.1 Dependence on the integration method 

It has been shown in Sect. 4 that the importance sampling 
integration method benefits from drawing samples in the 
extended space of the variables: 

{JVi > ... 1 JV lf |n > ...,rj V _i} with N = J2 N o (49) 

3 

instead of calculating the integral in the r's for each chan- 
nel {Nj} in Eq. l|2*ol with a fixed (say 1000) number of 
Monte-Carlo samples. One could try to apply the same 
idea to the Metropolis algorithm: instead of making a ran- 
dom walk in the space of channels {Nj} and evaluating the 
weight of the channel at each step by performing a Monte- 
Carlo integration of ^{Nj} m Eq. I|25|) . samples can be 
drawn in the extended space of the above variables and 
evaluating the integrand function <?" only one time at each 
step, thus saving CPU time. However, the relaxation time 
T re i will be different in these two cases and most likely 
longer in the extended space sampling case. This can be 
easily understood if the single \P evaluation is regarded 
as the extreme approximation of the integral in the r's in 
Eq. with N s = 1. This is shown indeed in fig. El 
the convergence to the equilibrium speeds up if the num- 
ber of samples used to estimate the integral in Eq. I(25|l 
increases. A serious drawback of the extended space sam- 
pling, i.e. with Ns = 1, is that the relaxation time may 
become so long that many bins of the multiplicity distribu- 
tion look like having reached their stable asymptotic value 
even when they actually still slowly drift towards it. We 
have checked this for the case shown in fig. I17l bv pushing 
the Metropolis random walk to 10000 steps, much beyond 
the scale of the step hystogram. It has been found that 
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Fig. 17. Step hystogram showing the convergence to equilib- 
rium of the probability Pe of a channel with 6 primary hadrons 
as a function of the step in a set of 500000 Metropolis random 
walks for different numbers of samples N$ in the Monte-Carlo 
integration of the phase space volume of each channel ^{Nj} 
in Eq. 125H . The horizontal solid line indicates the Pe value es- 
timated with the importance sampling method and the dashed 
band its relevant statistical uncertainty. Cluster mass is 8 GeV, 
energy density 0.4 GeV/fm 3 and charges are free. 

even at such large number of steps, the asymptotic value, 
calculated independently with the importance sampling 
method, is not attained, though the WOSSR test yields 
a positive response, under most circumstances. This find- 
ing indicates that the use of Metropolis algorithm requires 
more care than expected. At least a comparison between 
the apparent asymptotic stable values with different num- 
ber of samplings or a cross-check with independent calcu- 
lations is necessary. 




10 20 30 40 50 

Step 

Fig. 18. Step hystogram showing the convergence to equi- 
librium of the probability P2 of channels with 2 hadrons as 
a function of the step in a set of 100000 Metropolis random 
walks for the updating rule based on the MPD (bottom) and 
on a simpler method described in the text (top). Solid lines 
indicate the true value obtained by computing the phase space 
volume of all of the channels. The cluster mass is 2 GeV, the 
energy density 0.4 GeV/fm 3 and charges are free. The number 
of samples used in the numerical integration of Eq. 112"^ was 
Ns = 1000. For this comparison, resonances have been kept at 
a fixed mass and quantum statistics terms in ^?{jVj} have been 
neglected. 



just above or below. In the second case, a new particle 
is randomly chosen among all possible species. In the 
third randomly chosen particle of the current 

configuration is removed. 



7.2 Dependence on the updating rule 

In order to show the effectiveness of the updating rule 
based on randomly sampling the MPD, as discussed in 
Sect. 4, we have compared it with a simpler updating rule 
based on milder changes of the current configuration. This 
rule is as follows: 

1 . Three probabilities r/o , r/ + and r\- are chosen such that 

V0 + V++ V- = !■ 

2. A random extraction (0,+,—) is made according to 
the probability distribution defined by the 77's. 

3. Depending on whether the outcome is 0, +, — the over- 
all number of particles in the configuration is kept, is 
increased by 1 unit or decreased by 1 unit respectively. 
In the first randomly chosen particle in the cur- 
rent configuration is replaced with one having a mass 



For this updating rule, the proposal matrix T(m — > n) 
has been determined and optimal acceptance matrix set 
accordingly (see Eq. (|4*2"|l ). 

Indeed, this rule involves a slowing down of the con- 
vergence to equilibrium because the fraction of rejected 
transitions is much higher than in the MPD-based up- 
dating rule. This is apparent in fig. 1181 where the rele- 
vant step hystograms are shown for the multiplicity dis- 
tribution bin with 2 particles for a 2 GeV mass cluster. 
Whilst in the MPD-based updating rule the equilibrium 
is achieved within few tens of steps, in the above rule sta- 
bility is not achieved even after 3000 steps. Similar differ- 
ences are found for heavier clusters. Therefore, the MPD- 
based updating rule is much more efficient with regard to 
computing time. 
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Fig. 19. Relaxation times in a Metropolis random walk (see 
text for the definition) as a function of mass and charges of the 
cluster. Energy density is kept fixed at 0.4 GeV/fm 3 . Lines are 
drawn to guide the eye. 



7.3 Dependence on the cluster mass and charge 

We have studied the dependence of T re i on the cluster 
mass and charges, at constant energy density, by apply- 
ing the WOSSR tests to step hystograms of multiplicity 
distribution P n in different bins. For each mass and set of 
charges we have taken the highest T re i among the bins for 
which P n > 10~ 3 , as the relaxation time for that cluster. 
This study has been carried out for completely neutral and 
pp-like clusters; energy density has been kept constant at 
0.4 GeV/fm 3 . For each cluster 500000 Metropolis random 
walks (100000 for M = 10 GeV) have been performed up 
to 300 steps. The relaxation times are shown in fig. ^5] 

This just defined relaxation time shows an initial in- 
crease as a function of the cluster mass and drops there- 
after going from 8 to 10 GeV. For the present, we do not 
have a complete understanding of this behaviour. An in- 
crease as a function of the mass is expected because the 
number of channels sampled from the MPD whose weight 
is much larger than their correspoding microcanonical one 
increases owing to the persistence of the shape difference 
between multiplicity distributions in the canonical and mi- 
crocanonical ensemble (see figs. I12I13JI . as already men- 
tioned in Sect. 5. The drop, as well as the observed dif- 
ference between neutral and pp-like cluster, in the region 
8-10 GeV might be a genuine effect due to a local minimal 
distance between MPD and microcanonical distributions 
or an artefact of the cut on probability at 10 -3 . What it 
is important to remark is that the relaxation times are 
not larger than 0(100) in the region where microcanoni- 
cal calculations are necessary. In order to further improve 



Metropolis calculations, the same modification of the sam- 
pling distribution put forward at the end of Sect. 4 to re- 
duce statistical error in the importance sampling, could 
be carried over here so as to speed up the convergence to 
equilibrium. 



7.4 Comparison with the importance sampling method 

The performances of importance sampling method and 
Metropolis algorithm have been compared by studying 
the statistical error on the calculation of averages of sev- 
eral observables with the same computing resources. For a 
neutral cluster with 4 GeV mass and 0.4 GeV/fm 3 energy 
density, we have calculated the total average multiplicity 
and the primary multiplicity distribution 100 times, each 
time taking 100000 steps in both methods. As the sam- 
pled distribution at each step is the MPD in both cases, 
the used CPU time is approximately the same. The func- 
tion in Eq. (|25|l has been sampled one time per channel 
in both methods. 

As an example, we show in fig.lJOlthe results obtained 
for the estimate of the probability P5 of channels with 5 
primary hadrons. It can be seen that the Metropolis al- 
gorithm gives rise to a broader statistical distribution of 
the estimated values with respect to the importance sam- 
pling method. Moreover, the distribution is not gaussian 
and it is slightly asymmetric with also few cases of out- 
ranging estimates. On the other hand, the distributions 
for the total average multiplicity look quite similar in the 
two methods. It is worth pointing out that the a priori 
statistical error estimate obtained by using Eq. I|18|) for 
the importance sampling method is about 0.008, in good 
agreement with the found RMS of 0.0087 quoted in fig. 

M 

We have not investigated in much detail the sources 
of such a different statistical resolution, but we surmise 
that the ultimate reason of a larger error in Metropolis 
algorithm is the extra call to the random number gen- 
erator which, at each step, is possibly needed to accept 
or reject a proposed transition. Ultimately, this is an ad- 
ditional source of fluctuations in the Metropolis random 
walk which is absent in the importance sampling method 
where each extracted configuration is simply reweighted. 
Altogether, the importance sampling method seems to be 
better performing to calculate averages in the microcanon- 
ical ensemble. 



8 Conclusions 

We have calculated averages in the microcanonical ensem- 
ble of the ideal hadron-resonance gas including all light- 
flavoured known resonances up to a mass of about 1.8 
GeV. We have found that microcanonical average mul- 
tiplicities of different hadron species differ by less than 
10% from the corresponding canonical average (i.e. cal- 
culated by introducing a temperature) for clusters with 
relatively low mass, around 8 GeV, and energy density of 
0.4 GeV/fm 3 . This confirms and extends previous findings 
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Fig. 20. Statistical distribution of the Monte-Carlo estimates 
of the probability P5 of channels with 5 primary hadrons, for 
a neutral cluster with 4 GeV mass and 0.4 GeV/fm 3 energy 
density. 



|S] obtained with a restricted sample of hadron species. 
However, multiplicity distributions in the two ensembles 
show a clear difference in shape which seem to persist in 
the thermodynamical limit. Particularly, the dispersion of 
the total multiplicity distribution is, to a great accuracy, 
1/2 of the Poisson dispersion in the large volume limit; 
presently, we do not have any simple explanation of this 
fact and whether this is linked to our special choice of the 
energy density. 

A major point of this study is concerned with the nu- 
merical methods developed to calculate microcanonical 
averages: an importance sampling method, which has been 
proposed and used in this work for the first time; and a 
previously proposed Metropolis algorithm [7], whose per- 
formances have been improved by using as proposal ma- 
trix Poisson distributions with mean values set as those of 
the grand-canonical ensemble. They have been also em- 
ployed as sampled distributions in the importance sam- 
pling method. 

The Metropolis algorithm is capable to provide single 
samples of the microcanonical ensemble with unit weight 
and it is thus a suitable tool for event generation. Yet, 
it was proved to be less accurate than the importance 
sampling method for the calculation of averages for sin- 
gle clusters. Moreover, the Metropolis algorithm used for 
event generation requires much care, particularly a pre- 
liminary study of how fast the equilibrium condition is 
achieved. The convergence to equilibrium may depend on 
the observable to be analyzed and, more worrying, on the 
specific integration methods to evaluate the microcanoni- 
cal weights of the channels. 



However, the present study indicates that there is still 
room for a further improvement of the efficiency of both 
examined methods. An efficient way of calculating micro- 
canonical ensemble opens the way to test the statistical 
hadronization model at low energy and with respect to 
many more observables than those considered as yet. 
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Appendix 

A Calculation of the phase space integrals 

Here we summarize a method to calculate phase space 
integrals due to Hagedorn 9 . An analytical calculation of 
<P in Eq. H12|l can be carried out for two or three particles 
whilst for more particles the expressions get rapidly so 
complicated that a numerical computation is much more 
suitable. 

For two particles: 



$(M, mi, m 2 ) — 



47rp* eie 2 
T 2 ~M~ 



(A.l) 



where m,'s are the masses, ej = -\/p* 2 + mf the energies, 
M the cluster's mass, T the total available kinetic energy 
and: 



M 2 - 2{m\ + m|) + (m? - m 2 ) 2 



For N > 2 a function W of the momenta pi , . 
introduced: 



A' 



<2> = 



(4tt) 



t=l 

The function W reads 



. N N 

J A — 1 " A — 1 



(A.2) 

,Pjv is 

-,Pn) 
(A.3) 



,Pn) 



(47 



1 r r N 1 N 

* 3 (£piP*) ( A -4) 



where Pi = pj/pj, and can be calculated explicitely: 

W(pi,...,piv) = - ojV+1 lAT * 

x ^ (Ti . . .cr N {'^cr i p i ) N 3 

{<Ji...a N } 



(A.5) 



where a can be either +1 or — 1. This expansion involves 
a large number of terms even at relatively small N, so it is 
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more advantegeous under some circumstances, to calculate 
W differently. In fact, one can set in Eq. (|A.4|I : 



N x N 

shy. p« p ») = 7^j3 J d3u exp [ - 1 pipi ■ u ] ( A - 6 ) 



which leads to 0: 



JY 



I f°° 



sin(p. t u) 



(A.7) 



Setting now y = ufi, where fi is an arbitrary energy scale, 
to make integration variable adimcnsional: 



W(pi,...,pjv) 



1 



2ttV J 



11 (Pi/m 



2ttV J 



1 f l , (1 - x) N ^ A sm [f l- x 



dx 



r N-2 



n 



(A. 



where the last expression has been obtained by using the 
variable transformation y — x/(l — x). The last integral 
can be calculated numerically provided that the number 
of particles is at least of the order of 10, otherwise the in- 
tegrand function features strong oscillations which make 
most numerical integration methods failing. In this case, 
it is compelling to calculate W by means of its full expan- 
sion HA.5[1 . However, the CPU time needed to sum all of 
the terms rapidly grows with N, so that for N > 10 it is 
definitely preferrable to switch to the numerical computa- 
tion of the integral (|A.8J| (see also fig. [2J . The arbitrary 
scale [i can be set so as to obtain a good resolution in 
the numerical computation of W and we have found that 
l-i = l GeV is an appropriate choice for most cases at the 
energy densities we have been dealing with; if the number 
of momenta larger than 100 MeV is < 10, we set fi = 100 
MeV. 

In order to minimize the CPU time spent to calcu- 
late W by using its full expansion ljA.51) we have opti- 
mized the loop over all combinations {oi, . . . , cr/v} such 
that a j'Pj > by considering only those which are rel- 
evant. The momenta pi, ... , pjy are first sorted such that 
Pi > P2 > • • • > Ph and all possible iV-uples made of +1 
or — 1 are arranged in a multi- level tree structure shown 
in fig. |2] for the special case N = 4. If, for some iV-uple 
at some level Y^,f=i a jPj — 0' the same is true for all N- 
uples branching from it because of momentum sorting and 
of the tree structure, which is such that moving upwards 
one level a change in sign (+—►—) occurs. Note that 
replica A^-uples in fig. |^ are just meant to show the tree 
structure, though they are in fact not considered in the 
actual loop. 

The phase space integral in Eq. (|A.3(1 is now trans- 
formed by means of a sequence of changes of integration 
variable [7| pi — > ti — ► Sj — * Xi — > Z{ — ► r%. If t{ denote the 
particle kinetic energies and T the total available kinetic 




L ^\ +-+- 



Level Level 1 Level 2 Level 3 Level 4 

Fig. 21. Tree diagram to calculate W with for 4 parti- 

cles. 



energy: 



Pi = y/t l (t l + 2m i ) i = l,...,N 

U = Si — Si— % i = 1, . . . , N with sq = and s n = T 



Xi = ^ * = 1 JV- 1 



Xi+l 



i = 1, . . . , TV — 1 with xn = 1 



ri = z\ i = l,...,N-l 



(A.9) 



As a result of this sequence of transformations, the Dirac 
delta of energy conservation in Eq. (jA. 3f) is integrated 
away and the phase space integral reads: 

#= / d n ... f dr N - 1 r(n,...,r N - l ) (A.10) 
Jo Jo 



where: 



T(r u . . . ,rjv-i) 



(47r) iv T 



Nrp3-2N 



N 



,Pn) 



(N-iy. ±±< 

(A.ll) 

and pi , . . . , pn are to be calculated by going along the 
transformations (|A. 9|) for a given set (ri, . . . , rjy). In the 
form (|A.10|I ^ can be easily calculated by means of Monte 
Carlo integration. 



B Error estimates for importance sampling 

We have defined the best estimate of the mean value of the 
observable O in Eq. (|17|) . This estimator can be viewed as 
a random variable because the channels {TV,} are random 
variables whose distribution is 77 in Eq. (|20fl . Let us in- 
troduce the random variables f2, and fl taking values 
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i?jjVj}, 0({Nj}) and Il^.y respectively for a particular 
channel {Nj}. The following relations involving expecta- 
tion values hold: 



-u 



E 

{Nj} 



n 



{N 3 } 



n 



n 



{N S } 



Q 



n 



E " n = E O^Dtt^ 1 ^} = n (0) (B.l) 

where Eqs. © and (|1(J|> have been used. We can now write 
the estimator random variable (0) from Eq. (|17|) as: 



(0) 



^ fc=1 n fc 



(B.2) 



If the number of events M is large enough the numer- 
ator and denominator in the right hand side are gaus- 
sianly distributed according to the central limit theorem 
and the variance cr 2 ((0)) can be calculated by using the 
error propagation formula: 

i /v 2 M 

a 2 ((0» = <r 2 (N)^ + ^ 2 (D)-^ - 2cov(N, D)-^ (B.3) 

where N and D are the numerator and the denominator in 
Eq. (|B.2fl . and N and D their expectation values E^N) 
and Ejj(D) respectively. Because of Eq. <|B.1|1 : 

N = E/r(N) = MQ{0) D = Eii(D) = MO (B.4) 



Thus, by using Eq. I|B.4(I and the general variance defini- 
tions: 



«t(x) = E(x 2 ) - E(x) 2 
cov(x,y) = E(xy) - E(x)E(y) 



(B.5) 



we obtain: 

<7 2 (N) 

CT 2 (D) 

cov(N, D) 



1 

D 2 

TV 2 

N 



1 



AW 2 

(O) 2 
MQ 2 

(O) 
MQ 2 



Eii 

Ei7 
EiT 







n 2 

ft 2 

n 2 

ft 2 



I2 2 







n- 



(B.6) 



Therefore Eq. (|B.3I) can be rewritten as: 



^ 2 ((0)) = 



1 



MQ 2 

2 



E n [0 2 %) -{0) 2 Q 



+ (0) 
-2(0) 



-n 



n 2 

ft 2 







(0}Q 2 



(B.7) 



The estimator (II 71) is a biased one. This can be proved 
writing Eq. l|B~2)l as (0) = N/D and: 



E„l5 



D 

E/z(N) 



= E„(N)E„ ( i 



1 



cov(N, 1/D) 
o- 2 (D) 1 cov(N,D) 



Eii(D) E ff (D)3 



Ett(D) 



(B. 



using the definition of covariance and expanding the func- 
tion 1/D around its mean value. The bias B of the esti- 
mator (IB.2II can then be written as: 



(B.9) 



B = E n {{0))-{O) = E n (^)-(O) 



MJ? 2 



by using Eqs. (|BT8|) . flB^ , (|Tl4|) . 

Estimates of the bias (|B.9|) and the variance (|B.7|) can 
be obtained by replacing in Eq. 1|B.7(I the expectation val- 
ues with arithmetic means: 



1 M 

e ^mE 



(B.10) 



fc=i 



while (O) can be estimated through Eq. (|17() and the best 
estimate of Q is, according to Eq. (|B.1|I : 



q = e 



l M 

M ^ 

fc=i 



(fe) 



77 



(fe) 



(B.ll) 



The estimator (|17[) could be corrected for the bias by sub- 
tracting the estimate of i|B.9|l . However, the bias i|B.9(l is 
proportional to 1/M unlike the statistical error a^o) de- 
creasing like l/\/M and, thus, becomes negligible with 
respect to the statistical error for a large number of sam- 
plings. 



C Multiplicity distribution in the canonical 
ensemble 

The multi-species multiplicity distribution in the canon- 
ical ensemble with internal abelian charges Q, can be 
obtained from the generating function associated to the 
canonical partition function Z(Q). The probability of a 
single state in the canonical ensemble reads: 



1 



Z(Q) 

and the generating function 



e-^/^q^q (C.l) 



K 



which is essentially the Eq. (|18fl . 



G{\u ■ , Ak) = ^ E e~ E - JT ^Q II A f statc 

^ ' states j — 1 

(C.2) 
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where K is the number of hadron species and Nj s t a te the 
number of hadrons of species j in the state. The generat- 
ing function can be worked out by Fourier expanding the 
Kronecker delta in Eq. l|C.lj) 5]E]: 



G(Ai, . . . , Ai<) 

K 



1 1 



-4x7 f d M </>e lC ^exp 

*) M J--* 



Z(Q) {2 

E (2 ^ )3 1)F / d 3 P log(l±V^»-^^' 

(C.3) 



4=1 



where the upper sign applies to fermions, the lower to 
bosons. 

According to Eqs. (|CJ. 1|) and (|CJ.2|> . the probability of 
a given n-tuple of hadrons {Nj}, that is the multi-species 
multiplicity distribution, can be obtained by integrating 
the generating function over the unitary circle in the com- 
plex Xj planes: 



P({Nj}) 



K 



4=1 J A 3 



G(Xi,...,X K ) (C.4) 



Plugging Eq. (|C.3(I in the l|C.4|) and using the residual 
theorem, one obtains: 



P({N 3 }) 



1 1 



Z (2?r) 



M 



d M d> e lC ^ 



K 



j=l J CIA, 



E (Ti) nj+1 



- nj iqj -0 



A=0 

(C.5) 



where 2j( n ) is as in Eq. I|35(l . Now the exponential in Eq. 
IIC. 51) is expanded: 



exp 



n 



OO \ '*j 

(= F i)« J +ifi^!iil_2_ e -nmr<t> 



exp 



A 

n 4 



OO OO 



n e (*d 



^■^^A/ 



(rij + l)/i 4(raj) 4 



n,=l h„ =0 



OO OO 



zjj ,X/ 



li„ , =0 ra,- = l 



n ■ ' h„ . ! 



5^°=i nj/i ni = Nj. Therefore, the last series in Eq. I|C.6J) 
gives rise to: 



Njle-^ 1 *- 4, E (Tl) 



n 



'7i i 



(C.7) 



where {h n A indicates the set of partitions (in the multi- 
plicity representation) of the integers Nj, i.e. integers such 
that ^2^.-1 Tijh nj — Nj. The rij and hj indices actually 
run from 1 to Nj. Defining h nj — Hj and restoring 
the Kronecker delta, one can rewrite Eq. I|C.5(I by using 
the expression of derivatives in Eq. I|C.7(I as: 



P({Nj}) = 



K 



n e (td ^ n 

4=1 {fen,} "j 



6 4K) 



n ■ J hn ! 



<*q,e v -i ( c - 8 ) 



which is Eq. 



D Error estimates for Metropolis algorithm 

An estimator of the mean value of the observable O in a 
M -steps Metropolis random walk has been defined in Eq. 
(|44|l . As well as for the importance sampling method, this 
estimator can be viewed as a random variable (0) and so 
the values of the observables at each step 0( fe ). Therefore: 



a 2 «0)) = E«0) 2 )-(O) 2 = E 



M 



Thence: 



- {Of 
(D.l) 



pEE(0 (fe)2 ) 



M 2 



fe=i 



M 2 



Ee(o«o«) 



fe<i 



= -1e(0 2 ) + -^V E(OWOW) 

fc<i 

1 „ M M-fc 

= ^ e (° 2 ) + 4eE^) + <o> 



(D.2) 



k=l 1=1 



where A is the autocorrelation function. We can then 
write: 



h nK =o 

(C.6) 

Taking the derivatives of the latter expression with re- 
spect to the Aj's in Xj = according to Eq. 1C.5I) . one 
is left with non vanishing terms in the above equation if 



Ef =1 ow 



2' 



V M J 

E(0 2 ) M - 1 



M M-k 



M 



E E A ^ ( D - 3 ) 



fc=i i=i 
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If M 3> T auto , we can approximate the inner sum in Eq. 
l|D.3jl with the integral of the autocorrelation function R 
defined as: 

oo 

fl = 5>(Z) (D.4) 
and rewrite an approximate expression for Eq. (|D.3|) as: 



Eti ° (fe) 



a/ 



Therefore, Eq. IjD.ljl becomes: 



Am - <° 2) - ( ° >2+2Ji (d.6) 



which proves Eq. I)45[l . 
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